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ABSTRACT 


viii 


Integro-pertia.1 differential equations from a 
population "balance model for a continuous microbial pro- 
pagator in which rod- shaped organisms grow have been 
solved using the method of weighted residuals. Since 
the domain of interest is semi-infinite,^? linear combi- 
nation of Laguerre polynomials, which are complete in 
the space (0, co) has boon assumed as the trial solution. 
The orthogonality property of these polynomials has been 
used to advantage in securing an exact moment condition 
for tho population distribution. 


‘ The tridl solution is substituted in the pertinent 
differential equation and the unknown parameters in the 
assumed solution determined by orthogonalization of the 
equation residual with members of a family of functions 
complete in the space (0, co ) . Both the steady state and 
unsteady state solutions have been obtained using three 
different weighting functions and the results are in good 
agreement. ^ft is found that the .final solution is approached 
in both the cases with about eleven functions^ ^Sie solution 
to the steady state equation has been checked by two 
successive a cor oxime t ions. The batch propagator has been 



CHAPTER - I 


INTRODUCTION 

Population balance equations represent a number 
balance of particles on particulate systems in which the 
number of particles is continually subject to change by 
the disappearance of existing particles or appearance of 
new particles or both. In chemical engineering, they form 
a sound basis for the theoretical understanding of conti- 
nuous crystallization, grinding, microbial propagation, 
multi-phase suspensions and catalytic fluidized reactors- 
subjects hitherto handled with a good amount of empiricism. 

In all the above systems, each particle v_may be 
characterized by a specific set of properties, for example 
the size of the droplets in a dispersed phase reactor or of 
crystals in a continuous crystallizer or the physiological 
state of the micro-torganisms in a microbial culture. Eor a 
mathematical characterization of such systems, a population 
balance equation is found necessary. 

Whenever a system containing a large number of 
particles has to be studied., as in the Kinetic Theory of 
Gases, the overwhelming complexity of detail necessary to 
describe the state of the system and its interaction with 
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like the Boltzmann equation'., are usually partial differential 
and often in t egro— differential equations. The equations 
contain "Source" terms which account for the appearance of 
new particles and "Sink" terms which arise due to disappearance 
of existing particles. These "Source" and "Sink" terms dep en d 
on the system being described. 

1*1 The Need for Population Balance Equations : 

The applications of the principles of conservation of 
mass, energy and momentum in engineering contexts to homogeneous 
systems is familiar to the chemical engineer* In these cases 
the equations of change furnish an adequate mathematical 
characterization of the system being studied. However there 
are many problems of engineering importance involving hetero- 
geneous systems like two phase suspensions which cannot be 
handled within the framework of the conventional conservation 
equations alone; the proper formulation of the equations of 
change demands a population balance. Bor these systems, the 
distribution of entities about spatial (or mo* -likely property) 
coordinate axes is of dominant importance. Bor example 
continuous crystallizer design is based to a large extent on 
the size distribution of crystals in the vessel and catalytio 
fluidized reactor design on the activity and area distribution 
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give a better understanding of many select! cnal mechanism 
in growth. 

The general nature of population balance equations 
for systems in which an individual particle is characterized 
by a vector quantity has been given by various author s(2,3 } 9) . 

Hulb\jj?t and Katz(2) have formulated a class of problems 
in particle technology in terms of equations similar to those 
of classical statistical mechanics and shown how these 
equations can be tied in with the differential material and 
energy balances commonly used to describe the performance of 
pieces of chemical equipment. The main problems t reated 
are those of particle nucleation, growth and agglomeration. 

The assumption of simple forms of expression for growth and 
nucleation and absence of breakage or fission phenomena 
enabled the authors to simplify the original partial differential 
equations to ordinary differential equ ©tions in the moments 
of the distribution* Once the moments were known, the 
distribution could be fitted using the appnopriate orthogonal 
polynomials of the given space. The general equation given 
by them for a particulate system in which the individual 
particle is characterized by a vector x is 
3W_ 

— (x, t) + V. x W x ( x , t ) — h(x, t) 
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¥ = Non- Normalized distribution of particles in state 

between x and x + dx 
x = growth rate of particles 
Vis the gradient in x space 

h(x, t) is the net rate of appearance of particles in state x 

A. D. Randolph ( 3 ) has derived a general form of a 
population continuity equation. He has further transformed 
it into a more useful form applicable to the majority of 
engineering problems and focussed attention on the similarities 
with respect to the more common engineering conservation 
equations. 

1 . 2 Chemic a lEngineering Applications- : 

Randolph and Larson ( 4 ) have applied the idea of 
conservation of population for the analytical study of crystal 
size distribution in a continuous crystallizer. Equations 
relating population density to size, both in the steady and 
transient states have been presented and solved. Assumptions 
were made to simplify mathematically intractable cases. 

Crystal size distributions have been obtained for several 
modes of crystallizer operation - Arbitrary solids con- 
centration, seed crystal removal, product classification, 
and staged vessels. The effect of holding time and feed 
super saturation on crystal size has been investigated. 
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Behhken, Horowitz and Katz (5) Have given a mathematical 
formulation for determining particle size distributions in 
mixing vessels where each particle, while in a vessel grows 
or shrinks according to a prescribed differential law. The 
application studied is a catalytic fluidized reaction system 
where the catalyst particles are recycled between two zones, 
one in which they are loaded with reagent „ and one in which 
they are stripped of product. Steady state size distributions 
have been developed along with expressions for average loading 
and stripping rates. Procedures are given for following 
transients. 

Fredrickson, Tsuchiya and Co- workers (6, 7,8) have 
introduced the idea of Conservation of population with 
respect to age and size distributions in microbial cultures. 
Fredrickson ct.al have also provided more general equations 
for microbial population - by attributing a "vector state" to 
characterize an individual cell (9) . In this"vector state 
description" of the microbial cell, the physiological state 
of the cell is specified by its biochemical composition, 
assuming that the intr a- cellular structure is of a relatively 
simple nature. Growth and fission are described by a Markov 
process, that is these depend only on the cell's current 
physiological state and not on how that state had been 
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derived, as well as an equation of change 'for the state 
of the cellular environment. These equations allow us 
to predict the statistical and dynamical behaviour of a 
cell population from information obtained by analysis of 
cellular and sub-cellular structure and function. 

1 • 5 Solution of Population Balance Equations ’ 

Graudin and Meloy(ll) have derived a comminution 
distribution equation considering both single and repeated 
fracture. The resulting integro-differential equation 
has been solved using a finite difference approximation. 

K.J. Reid ( 1 2 ) has presented a practical approximation 
and developed a. solution for the fundamental integro- 
dificrential equation for batch grinding. The solution 
has been written in terms of two experimentally determined 
parameters of the system - rates of breakage and a breakage 
function. 

Valentas 'md imundson ( 1 3 ) have developed a mathe- 
matical model which relates the breakage of droplets in a 
two-phase system to the distribution of droplet sizes. 

The conservation equation leads to a Volterra Integral 
equation in the steady state. This involves the influent 
distribution function, the vessel distribution function 
and a Kernel describing the breakage mechanism. This has 
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feed distribution, perfect mixi and binary breakage 
with formation of equal sized daughter droplets, i 
numerical solution has been developed for the more compli- 
cated case where a normal distribution about the mean 
of the sizes of the daughter droplets formed from breakage 
of larger droplets is assumed. The integral term is 
approximated by an integration formula and the solution 
of the equation obtained by collocation Method. Solutions 
have been also obtained for the case where the breakage 
frequency depends on size and for limited breakage where 
droplets below a certain size are stable and do not break. 

In another paper (14) where coalescence is also considered, 
an expression for the coalescence frequency distribution 
has been developed by an analogy to collision in a second 
order reaction. The steady state equations have beer- 
solved numerically with an iterative procedure for the 
final system of nonlinear eciuations. The transient equations 
for a stirred-tank vessel and a hatch vessel have been 
solved by approximating the breakage and coalescence 
integrals by integration formu 3 .ee and solving the resulting 
system of differential equations by a predictor-corrector 
method. 

Ramkrishna (lO) ot.al have solved the population 
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of successive approximations on a computer. 

Eakman ( 7 ) has solved the population balance for 
a continuous microbial propagator based on a mass model 
for the case where the cells divide into two identical 
daughter cells. Solutions were obtained by a predictor- 
corrector method for spherical and rod-shaped cells. 
Exponential batch propagation is considered on similar 
lines. 

1 .4 Scope of the Present Work : 

Population balance equations become integro- 
differential or integral equations when "new" particles 
are formed from existing particles or due to interaction 
between existing particles. .Additional complication is 
introduced when there is growth, as in a. biological culture 
with which the present work is concerned. Growth terms 
can arise when any property associated with the individual 
particle changes with time. Thus in dispersed phase 
systems if one chooses to describe the mean concentration 
of a transfer species as a property of the individual 
droplet, then this would change with time as a result of mass 
transfer . 


Experiments may dictate a maximum size of particle 
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organisms for example, the form of the growth rate expression 
does not permit one to fix a maximum size of cell. Hence an 
infinite domain may arise. 

Solutions to population balance equations of the 
above type are non-existent in the engineering literature. 

It is intended to solve the population balance 
equation for a continuous microbial propagator in which 
rod-shaped organisms grow by a classical approximation scheme, 
the method of weighted residuals. 
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CHAPTER - II 

POPULiPIOE BALANCE BQU-APTOHS FOR A MICRO- 
BliL PROP^GiTOR 


The system with which the present work will he 
concerned is a perfectly mixed vessel containing microbial 
cells and substrates. i feed stream containing fresh 
substrate but no cells is introduced into the vessel, while 
an effluent stream of the same volumetric flow rate removes 
the cells, substrates and products from the vessel. A 
constant culture volume is thus maintained within the vessel. 
Phis system will in future be referred to as a "continuous 
microbial propagator". 

Phe cells in the propagator grow by consuming substrates 
from the environment and multiply in number by a random binary 
fission process. Growth and increase of population bring 
about changes in the cellular environment both by uptake and 
release of chemical products by the cells. Phe nature of the 
cellular environment in turn affects the rate of growth and 
multiplication of the population. Phe mathematical model 
adopted for the description of the behaviour of the microbial 
population in the propagator will be the one developed by 
Enkmarif 7. B) . 
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The model is segregated in the sense that cells are 
distinguished as individuals* Further, it is structured, 
so that differences between cells are recognized. The 
interaction of the population with the environment will be 
explicitly taken into account. Bakman ( 7 ) considered both 
cell mass and cell age and concluded that cell mass is a 
better index of the state of the cell. Hence the mass model 
will be used throughout. 


2.1 The Population Balance Equation for a Continuous Microbial 


Define 

W(m, t) dm = the number of cells per unit volume (population 
density) with mass between m and m + dm 

P = ' tlle probability a cell of mass m at time t will 

divide in time t to t + dt. 

p(m, m')dm = the probability a daughter cell formed from a 
mother cell of mass m’ has mass between m and 


d + dm. 

f 

G s j c )dt = the probability a cell of mass m at time t 
will die in time t to t + dt. 


r C slP 


the rate of mass increase (growth rate) of a 
cell of mass m in environment with nutrient 


( n 
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Ihe partial integro differential equation describing 
the distribution of masses at any time is obtained by making 
a number balance on cells of a given mass, similar to the 
technique one adopts in deriving the conventional equations of 
change, or that given by Randolph ( 3 ). This gives 


^ W (m, t) 
> t 


+ 



[r(m, C gfc ) W(m, tf} 

(m*, C gk ) p(m, m l ) W(m l , 


t) dm’ 


i + r (m > <w + 


W(m, t) 

( 2 . 1 ) 


The first term on the left hand side represents the accumu- 
lation term. The second term is that due to growth in and out 
of the given mass range. The first term on the right hand 
side refers to the number of cells formed by binary fission 
of larger colls to produce daughter cells in the given mass 
range* The second term includes the loss in number by washout, 
fission and death respectively? 0 is the holding time, the 
transition probability for washout in a perfectly mixed vessel 
being At/0 . 

Boundary conditions on the solution are obtained by 
considering an independent balance on population density in 
the propagator 

co r 1 

"i *vr / "vr / s*"***. 
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The number density N(t) is related to the distribution 

Wfa, t) as its zeroth moment. Thus 

no - 

IT(t) = J W(m, t) dm ...(2.3) 

Q 

Integrating equation (2.1 ) over all m and comparing the 
resulting equation with (2.2) one obtains 

r(co , C gk ) W(co , t) = 0 ...(2.4) 

which is the regularity condition stated by Behnken et.al.(5). 


Since the behaviour of a cell depends upon the substrate 
concentration in the environment, a substrate balance is made 
on the propagator. Bor an arbitrary ith substrate a mass 
balance gives 


dC 


00 


si 


dt 


3 (C si ~ C si ) “ jA (m) r ' (m ’ C sk } W(m ’ dm 


OO 


r n (m) W(m, t) dm 


...(2.5) 


Where C° . is the concentration of the ith substrate in the 

o -L. 

feed stream; Z 3 ^ represents fraction of ith substrate in the 
total amount of mass taken up the cell; r’(m, C gk ) is the 
rate of mass tip take of a cell. The net rate of mass increase 
for the cell r(m, C gk ) is the difference between the rate of 
mass uptake r’(m, C gk ) and the rate of mass release by the 
cell r"(m) 
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/ **f i is the fraction of ith product in the total mass released 
by the cell. 

Equations (2.1) and (2.5) subject to the boundary 
condition ( 2 . 4 ) and appropriate initial conditions thus form 
the mathematical problem which is to be solved. 

2 • 2 Growth of a Si ngle R o d- Shaped Cel l- 

Von Bertalanffy (15,16) reasoned that since all mass 
added to a cell must cross the cell surface, the rate of 
uptake should be proportional to the surface area of the cell. 
Further the rate of release should be proportional to cell 
mass since materials released are presumably products of 
catabolic reactions. It was assumed that the rates of these 
reactions were proportional to cell mass. Thus 


= r(m, C sk ) = 0'S - ...(2.6) 

By knowing how the cell geometry changes with growth 
for a rod shaped cell, an expression may be written for the 
growth rate ( 7 ). 

J A 

...(2.7) 


dm 

dt 




c ' 3 

j* is the volume average density of a single cell and R the 


cylindrical radius. It is assumed that the flux term^S is 


of Michael is-Menten form 


(c ) 


A 


s 


...( 2 . 8 ) 
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where is the concentration of the limiting substrate in 
the environment. It is assumed henceforth that there is a 
single limiting substrate. The second term in (2.7) is the 
contribution of the spherical end caps of the cylindrical 
cell and may be usually neglected. 

Integration of (2.7) for constant environment provides 
after neglecting the end ca,p contribution, 


m = m. 




...(2.9) 


Equation (2.9) shows that rod-shaped cells increase in size 
in an exponential manner with cell age and can in a mathe- 
matical sense reach arbitrarily large sizes. Thus the domain 
of interest should be (0, oo ) for rod-shaped cells. This is 
not so for spherical cells (7) where there is a limit to the 
maximum size attainable. 


2.3 Single Cell Division ". 

Schaechtor et.al(l7) have found that the size of a 
cell at division has a smaller coefficient of variation than 
cell age at division which indicates that cell size and hence 
cell mass is a better predictor of cell division than cell 
age. Cell division being a random process, it is found that 
the masses of cells at division are found to be distributed 
around a mean. A less realistic assumption would be that 
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assumed in this work to he of a Gaussian type 


h’(m) = - — ?~ e ...(2.10) 

€ Ja £ 1 + erf (~jr) ] 

By Bayes theorem of probability, one can relate the probability 
that a cell of mass m will divide in time t to t + dt, p(m,C g )< 
to the probability that a dividing cell has mass m, h*(m)« 


P (m, 0 ^) at = 


m - mc \ 2 
2 e ^ e [ 

erfc( 


~y (m , c sk) 

m - me \ 
£ ' 


dt 

...( 2 . 11 ) 


using dm = r(m, C dt 

£ is the standard deviation of the distribution. Since 
r(m, C g ^) depends on cell geometry and time, {"""(m, C^) will 
also depend in general on geometry and time. 

In the cell division process, a large number of repli- 
cated entities like nucleic acids are allotted to each daughter 
cell. If this allocation is random, then the mass of the 
newly formed cell will be random. In the present work the 
distribution of the cell masses at division will be chosen 
arbitrarily. One can choose any distribution satisfying the j 
symmetry condition p(m, m’) = p(m T - m, m’) and the normal!-* 
zation condition 

,-W 

l p(m, m’) dm = 1 


...( 2 . 13 ) 
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ftince a cell a.t division cannot have a mass more than the 
parent cell. It is assumed that the distribution of the 
daughter cell masses does not depend on the cellular environ- 
ment. The distribution of daughter cell masses that will be 
used here is 

p(m, m’) = ...(2.14) 

m ’ 5 

which satisfies (2.12) and (2.13) 


2.4 The Biomass Bal a nce f or a Continuous Propagator : 

On neglecting the spherical end cap contribution to 
the growth of rod-like cells, equation (2.7 ) can be written as 
2 0(0 


dm 


XT 


- A 




m 


...(2.15) 


The first moment of (2.1) is 
co 
dC 


dt 


r(m, C gk ) W(m, t)dm - — 

oo 


m 


(m, C sk ) ¥(m, t) dm 


..*(2.16) 


where the first moment 

sGO 


A, it) 


m W(m, t) dm 


. i . (2. 17) 


j 

o . 

is the biomass concentration. Substitution of (2.15) into 

(2. 16) and (2.5) and integration gives for the rate of change 
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dC _ 
dt 




A 


I +~c“ 


s 



. .. (2.18) 


and for rate of change of substrate concentration 

C C 

+V) + ^Ac C 

si 

...(2.19) 

Equation (2.18) is a condition for the first moment of W(m, t) 
that can he determined once the cell parameters are known, by 
simultaneous solution of (2.18) and (2.19). It must he noted 
that the form of (2.19) leads to biomass and substrate balances 
(2.18) and (2.19) which do not involve W(m, t) . This is not 
true in general where (2.15) could have a more , complicated form ; 
as for spherical cells where biomass and substrate concentration 
cannot be obtained without first solving for the cell mass 
distribution, alternately if the kinetics of growth is known 
in terms of distributed model parameters like biomass and 
substrate concentrations then it is possible to have the first 
moment and substrate balance eauations uncoupled from the 

i 

population balance equation. 

Equation (2.15) can thus be cast in a form . . 

r(m, C ) = f"( t) m where f(t) contains substrate concentration 
s 

terms which are determined from equations (2.18) and (2.19). 


dC s 

dt 


1 

e 


/ ° 
<°a ' 


°s> * 


A 2 

R P 


/r~ 

V s 


Equation 2.1 becomes 



19 


j^_W(m, t) 
t 



t ' ) dm f 
W(m, t) 


...( 2 . 20 ) 


Equation (2.24) has to Toe solved for the mass distribution 
W(m, t) , together with eouations (2.22) and (2.23) which 
determine f(t) and substrate concentration C g , subject to 
the initial conditions 


and 


W(m, o) = g(m) , an initial mass distribution 

f-Oo 

C(o) = I mg(m)dm 

C (o) = an arbitrary initial substrate concentration. 

o 


2.5 Various Modifications of the Problem ; 

Eirst, the steadv p+ate problem will be considered. 
Equation (2.20) is written for the steady state 

~ c /v/ 

|jr(m) W(m)J t= 2 Jp(m’, C g ) p(m, m T ) W(m’)dm l 


-& 




W(m) 


,..(2.2l) 


assuming that there is no death. mf(t) has been replaced 
by r(m) where r(m) is the steady state rate of mass increase 


for a single cell. 

r 


■(m) = 


2 

TT7T 


m 


rs/ 


- a! 


m 


.. .(2.22) 
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C is the steady state limiting substrate concentration, 
s 

Substitution of (2.22) in eauation (2.16) for the steady state 
yields 


( 2 

R f 


0 


0 

m s 


K + C 
8 £ 


- /o - 1 


© 


and (2.22) becomes for the steady state 


•(*) = 


m 

e 


. .. (2.23) 


where (? is the holding time of the reactor, it steady state 
simultaneous solution of equations (2.18) and (2.19) with 
their left hand sides equated to zero gives the steady state 
substrate concentration 


^ _ (l +Q/ C c) 

2 Q 

-Pr~ « 

and the steady state biomass concentration 


...(2,24) 


~ c: - 
c = - 


e 


— 


(2 1 fin. ± 


f ' K + 'C 
s s 


^Ac] 


...(2.25) 


It is assumed that the limiting substrate is not given off 


by the cell, 


1 


= 0 . 


Several interesting variations of ( 2 . 2 1 ) arise depending 
he 

p(m, m') 


upon the nature of the probability functions f '1 (m, G g ) and 


'* , The simplest model would consider division of cells 
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cells dividing into daughter cells of identical mass (m_/ 2 ) 

c 

In such a case the probability that given a cell has divided, 
mass is m is h^m) = £ (m - m c ) a Dirac delta function. 

(m, C s ) = o (m - m c ) r(m, C g ) 

The probability function p(m, m f ) = £ (m - ~) equation ( 2 . 21 ) 
becomes 


dW 

dm 


r(m, C s ) 


JL + dr 

0 dm 


Qw(a) 


. .. (2.26) 


Dhe solution (7) "io (2.26) using the moment condition 

pQ 

C - j m W(m)dm 

a 

to evaluate the integration constant, is 

~ m c 


W(m) = 


— 2 ~ { ni ( 


m In 2 
W(m) =0 m >• m Q 

m < i m c 

This model is a gross oversimplification. 


(2.27) 


II Cells divide on reaching a critical mass m , that 
is |“'(m, C g ) is still a. delta function, but the daughter cell 
masses are distributed according to 

P (m, m') = lojai i al-jah: ...-(2.14) 

m’ 5 

Equation ( 2 .' 2 l) becomes in this case. 
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OO 
f 


d_ 

dm 


[w(m) r(m, C s )^ = 2 |r(». 0 B ) <f(»' - ^ *<»’> P<»-“’ )aa ’ 


W(m) 

e 


or d_ j^ r ( m ) w(m)] = 2 W(m c ) p(m, m c )“- " -*^ 2 * 28) 

Using (2-22) for r(m) , an analytical solution can fee written 
down for (2*28) 

/yvv 

. . .(2.29) 


m^W(m) = m c 2 W(m c ) ^ m p(m ? m c )dm 
Let Kj = 2m c W(m c ) 


7Q m 4 m 
W(m) = Kj ^ \ — — - 

m / L 4 


5 6 

2 m_ nr , m_ 

5 c + 6 


on using (2.14) for p(m, m c ) in (2.29) and integration. 
The integration constant is evaluated using 


OfVV 


c = 


■c, 

m W(m) dm 


O ^ 

2 C 

K ! = TT5T C 


and 

W(m) 


60 0 

37 m 


[• 


15 m 2 m„ - 24 m m 3 + lOm 4 ^ 


for m <( m c 
for m *> m 


. ..(2.30) 


W(m) = 0 

This solution represents a tetter model than X-' though still 
an oversimplification. It is nevertheless useful in serving 
as a comparison with the more complicated models to he dis- 
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approximation to the actual solution of equation ( 2 , 21 ) in 
the regions of small mass where both p(m, m l ) and p (m, C g ) 
would be expected to be very small. 

Ill The masses of cells at division are distributed 
according to equation (2.10), a Gaussian distribution. Here 


p u, \^) = 


2 e 




m - m 


c \ 2 




r(m, 


m - m 

®2?fc ( — g— ~ ) 


... ( 2 . 11 ) 


The daughter cells are identical in mass that is 
p(m, m T ) = (m - |L) , Equation (2.21) becomes 

W ^) = 4 p (2m) W(2m) ~ (j. + f (m) 3 W(m) 

... (2.31) 

a functional differential equation, solutions to this 
equation have been presented by Eakman (7,8), by using a 
predictor-corrector method. The solutions to (2.31 ) by the 
method of weighted residuals followed by a successive 
approximation procedure will be giveh in Chapter lVi This 
model -is more realistic than (l) or (2). 


MV It is the object of this work to solve the most 
general of the above models in which distribution of cell 
division masses as well as distribution of daughter cell 
masses are considered and described by equations (2.-10) and 

l O HP A 1 -M Jm. mm. A 4- •£ A ** * . „ _ ^ . 1 * _ 
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by the method of weighted residuals for the steady state 
equation ( 2 • 2 1 ) as well as the transient equation (2.20). 


2*6 Batch Propagation ; 

The equations of the cell model for batch propagation 
are obtained from (2. 20), (2.18) and (2.19) by setting the 
dilution rate 1 /^ = 0* The partial integro differential 
equation describing cell mass distribution is described by 


c© 


■& w , ^ p 


jjr(m, C s ) W(m, t)j = 2 i p (m { C s ) p(m, m l ) W(mpt)dm' 


KYX 


- p(m, C s ) W(m, t) 


...( 2 . 32 ) 


and the substrate and biomass concentrations by 

oO 

dC 


s 


dt 


r- , f — 

= |oL(m) r"(m) - (m) r ! (m,C g ) 

at 


o 


W(m, t) dm 

...(2.33) 


dC 

cTT 


oO 

r 


r(m, C g ) ¥(m, t) dm ...(2.34) 

assuming no death. (2.32) is to be solved subject to (2.33) 
and (2.34) and the regularity condition 

r(co , C g ) ¥( co , t) = 0 


2.7 Exponential Phase Propagation : 


Consider a batch culture of a single species propagating 
in the exponential phase. If it is assumed that the substrate 
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concentration is sufficiently high, the changes produced by 
growth have no noticeable effect on the substrate concentration 
and hence no effect on the rate of population increase. If 
the culture has been growing long enough so that the cell 
mass distribution has become independent of time, the cell 
mass distribution takes the form 

W(m, t) = exp ( 9 t) g(m) ...(2.35) 

where N 0 is the population density at some arbitrary zero 
in time in the exponential phase. 


Substitution of (2.35) into the number balance (2.32) 
gives the integro differential equation for the asymptotic 
batch distribution 


dg* = 2 

dm r(m, C ) 

o 



1 

r(m, TT) 


C g ) g(m') p(m, m') 



C s ) + 9 + 



g(m) 

...(2.36) 


It is seen that equation (2.36) is exactly similar to the 

steady state continuous propagator population balance (2.21) 

except for the appearance of the term )) instead of (1 /q ), 

substitution of (2.35) in the total number balance equation 

(2.2) applied to a batch system provides a relation for , 

the specific rate of population increase. 
oO 

9 = ( fl (m» Cj g(m) dm 


...(2.37) 
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To calculate -V the distribution has to be known* Therefore 

(2.36) and (2.37) have to "be solved simultaneously for g(m) 

along with the normalization condition 
ctO 

g^m)dm = 1 ...(2.38) 

Once again the various models discussed for the continuous 
propagator would be applicable here also, subject to the 
additional condition (2.37). 




CHAPTER - III 


THE METHOD OP WEIGHTED RESIDUALS (M¥R) 

The method of weighted residuals is one of the powerful 
tools used by the engineer to find approximate solutions to 
the equations of change of distributed systems. The approxi- 
mate solutions in many cases provide a reasonable first 
guess that can be used in a scheme of successive approxi- 
mations to obtain improved solutions. The analytical form 
of the solution obtained by MWR is very useful if there are 
many quantities of interest to be derived from the solution 
to the problem. The method is applicable to linear as well 
as nonlinear problems and requires much less computation time 
to generate the solution than finite difference techniques. 

Collatz ( 1 8) has covered the method in detail, with 
many examples, applied to ordinary and partial differential 
equations. He calls them "error distribution" principles 
in accordance with the notion of distributing the error as 
uniformly as possible throughout the domain of the solution. 
Crandall(l9) and Kantorovich and Krylov (20) have extensively 
applied the method to problems in the mechanics of solids. 
Einlayson and Scriven ( 21 ) have given an excellent review of 
the method and its application along with an exhaustive 
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some aspects of the method to other methods like separation 
of variables and variational methods has been outlined. 

3. 1 Basic Method for Initial Value Problems - 

Let the mathematical model of a distributed system be 
represented by a system of differential or integro-differential 
equations. The boundary conditions represent the interaction of 
the system with its environment and the initial conditions 
represent some base state of interest. The general procedure 
in MWR is to assume a trial solution containing an arbitrary 
number of known functions of position and unknown functions of 
time. The undetermined parameters are found by substituting 
the trial solution into the differential equation and requir- 
ing that the error involved be distributed over the solution 
domain in some specified manner. 

Consider the differential equation for W(x, t) 

= IT(W) ~x in V, t 7" 0 ...(3.1) 

where N( . ) denotes a general differential or differential 
integral operator involving spatial derivatives and integrals 
of W. V is the domain of interest with boundary S and. t 
represents time, let the initial and boundary conditions be 

W(x, o) = W (x) x in V 

W(x, t) = f (x, t) x on S 

D 


.. .(3.2) 
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.Assume a trial solution of the form 

N 

w*(x, t) = 0 g (x, t) + ^ChCt) 0j_(x, t) ...( 3 . 3 ) 

1=1 

where the approximating functions 0 i are known and satisfy 
the boundary conditions, 

^s ^ f s ’ 0i = 0 x on S ...(3.4) 

In this case, W* satisfies the boundary conditions whatever 
Ch(t) may be. Ihe differential equation residual and the 
initial residual 

R(W*) = N(W*) - -|t* ...(3.5) 

^ N 

R 0 (W*) = W Q (x) - 0 g (x, 0 ) - 2.0 i (o) 0.(x, o) ...(3.6) 

-tsi 

indicate the extent to which the trial solution satisfies the 
differential equation. For an exact solution, the residuals 
vanish identically. As the number of approximating functions, 

N is increased, the residual could be expected to become smaller. 
Ihe residual is made zero in an approximate sense by equating 
its weighted average over the whole domain to zero. 

<^r(w =0 ...( 3 . 7 ) 

<^R 0 (W*)/f\j^ =0 3 = 1, 2, N 

where 

O v > s S 

V 


uvdv 


.. .( 3 . 8 ) 
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chosen weighting function. Equations (3-7) become a set of 
ordinary first order differential equations in the N unknowns 
C^(t) with the second of the equations providing the initial 
conditions C^(o) . 

W 

Eor a linear problem = l(¥) equations (3.7) become 


N , n N 

<Tr ^i> = ^ ...(5.9) 


i = 1 

i = 1 



I( V> 


or, in matrix form 



A f! = io + i 


...(3.10) 


Subject to C(o) = C Q 

The solutions to this system provide the approximate solution 
(3.3) of the problem (3.1 ) .Improved approximations are obtained 
by increasing N and their convergence should give an indication 
of the approach of the approximate solution to the exact solution. 
In the above formulation, the boundary conditions were 
satisfied exactly by the trial solution and the differential 
equation approximately. This is called the "interior method 1 '. 

We speak of a "boundary method" when the trial solution 
satisfies the differential equation exactly but not the boundary 
conditions. When neither the differential equation nor the 
boundary condition is satisfied exactly, we have a "mixed 



31 


is made depending on the control required over the relative 
accuracy of the solution at various points in the region. 

A mixed method would allow flexibility in the choice of trial 
functions. 

3.2 Choice of Approximating Functions 0. 

- 4 — - - - -- jL 

The choice of approximating functions is primarily 

dictated by the domain of interest. The functions must be 

linearly independent and differentiable to the extent that 

all terms in the differential equation and boundary conditions 

can be obtained. If the system of functions is complete, 

that is any piecewise continuous function in the domain can be 

approximated in the mean by an infinite linear combination of 

the basic functions, then convergence of the method may be 

expected. Mathematical proofs of convergence (20, 24» 25) are 

available for some problems and allow the determination of 

various sets of functions 0.* that are complete for the given 

J 

problem. If the 0^'s are chosen as the first n members of 
such a set, the approximate solution must approach the exact 
solution as n tends to infinity. Kantorovich and Krylov (20) 
state that the sufficient condition for convergence to the true 
solution as we increase the number of trial functions is the 
completeness of the system of families of trial functions. 

When the region of interest is finite, the 0^ are 
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p 

members of a family of elementary functions like 1 , x, x , .... 
or Sin x, Sin 2x, Sin 3x, ...etc. The choice between families 
rests on the nature of the problem at hand; if one were con- 
cerned with a vibration problem, it would be natural to choose 
the latter of the above mentioned families. Families with 
polynomials as members are popularly used as trial solutions. 

In these cases by using approximating functions which satisfy 
the boundary conditions and have the proper symmetry properties, 
accurate results can often bo obtained with a very small number 
of functions. Snyder, Spriggs and Stewart (22) have chosen 

p p p n p 

approximating functions of the form (C - x ), (C - x ) , ... 
for obtaining velocity profiles for the steady state flow 
between two inclined plane walls where a boundary condition 
must be satisfied at x = + C.- Gne finds this technique used 
very often in the Pohlhausen method of boundary-layer theory. 
Periodic functions have been used by Snyder and Stewart (23) 
to obtain velocity and pressure profiles for Newtonian creeping 
flow in a regular packed bed of spheres. 

Very few problems have been reported in the engineering 

literature using MWR for an infinite or semi- infinite domain. 

It would seem natural to use elementary functions belonging to 

the family e“ x x 11 or e'" nx as approximating functions for the 

2 2 

semi- infinite domain and e~ x x 11 or e~ nx for the infinite 
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Orthogonal polynomials as aporoximating functions are 
very convenient to use in many instances. For example, when 
one is solving for a distribution function as in the present 
work, expression of the trial solution in terms of orthogonal 
polynomials permits ono to write the leading moments of the 
distribution straightaway in terms of the coefficients of the 
trial solution expansion. If a distribution W(m) can be 
expressed as 

ET 

W(m) = C p(m) 0 (m) where 0 (m) are orthonormal 

ii n n 

n=0 

polynomials with respect to a weight* function pCm) , then the 
inner products with respect to p(m) , over the domain of the 
distribution, 

<W(m) , 0 o (m)^> = C Q 

<^W(m) ? 0 1 (m)/ > = C 1 ...(3.12) 

< \W(m) ^ 0jy(m)/ > = C]\j 

0 n (m) is an nth order polynomial and the system (3.12) can be 
written using the orthonormal property of 0^ as 

K oo*o = °o 

K ioA/ K n.Ai = °i 
N 

S K N;j/t) = 
j=o 


...( 3 . 13 ) 
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where K. . is the coefficient of m 3 in 0.(m). System (3.13) 
is a closed system and can he solved recursively for the 
moments JA- 0 , in terms of coefficients of the trial 
solution C Q , C^ , Cg, . ...C^. Alternatively, using (3.13) 
the trial solution (3.3) may he written in terms of moments 
and condition (3.7) would generate equations in the moments of 
the distribution. The utility of MWR seen in this regard is of 
immense significance to problems involving distribution functions. 

Orthogonal functions suitable for the finite domain 
could for example be Jacobi or Legendre polynomials, for semi 
infinite and infinite domains, Laguerre (See Appendix A) and 
Hermite polynomials respectively are found suitable. The 
subject of orthogonal polynomials is extensively treated by 
S*ego (27). 

3.3 Choice of Weighting functions h/''. 

- , ' J- 

The choice of weighting functions , corresponds to 
various criteria in MWR. 

(i) Collocation Method * The residual R(W*) in equation (3.5) 

is made to vanish at R collocation points. The weighting 

function is a Dirac delta function, *h . = & (x- - x) 4 

J d 

j = 1 to IT. The collocation points are usually chosen at 
regular intervals, 

A modified collocation method has been suggested by 
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function as well as the collocation points are determined 
by the orthogonality of the residual to the trial functions 

themselves. The method is applicable easily when the resi- 

« 

dual can he expressed as the product of a polynomial times 
a function. The advantage of the method is that in addition 
to averaging the residual to zero over the whole domain, it 
is made zero at the collocation points as well. 

(ii) Subdomain Method : The domain of interest is divided into 
N sub domains and the weighting functions are 

=1 x in V. 

J *1 

=0 i not in V. 

J 

j = 1 to N 

is N is increased, the sub domains become smaller and 
smaller and the residual is approximately zero everywhere. 


(iii) Least Squares Method ; By choosing the weighting function 

j = 1 to IT 


^3 ■ 


~i> R(w*) 


the mean square of the residual 


T 


(R(W*)) 2 dV 


V 


is minimized with respect to the Ch. 

(iv) G alerlrin's Method (20) : The weighting functions '/'h 
are the same as the approximating functions 0^. If the 
approximating functions are members of a complete set of 
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the residual orthogonal to members of the complete set. A 
complete set of functions has the fundamental property that a 
piecewise continuous function can be orthogonal to each and 
every member of the set only if the function is identically 
■zero. If the approximating function and the differential 
equation are such that the residual is continuous, it can 
vanish only if it is orthogonal to each member of a complete 
system of functions* In practice, the residual is made ortho- 
gonal only to a finite number of the members of a complete set. 

(v) Method of Moments : Galerkin's method is really a particular 
case of the more general method of moments in which the residual 
can be orthogonalized with any complete set of functions and 
not necessarily the trial functions themselves. G-alerkin's 
method was singled out only because it could be related to the 
Rayleigh- Hits method (20) used in solving variational problems. 

The method of moments has been used throughout the 
present work. It was felt that a comparison of the solutions 
of the problem with different weighting functions would give 
a good indication of the error involved in the approximation. 
(See Chapter IV) . 

The special case of the method of moments with weighting 
functions identically equal to 1 corresponds to the Rohlhausen 
method of boundary layer theory. 
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3.4 KWR vs Other Methods for the Given Problem: 

For the problem at hand MWR seems to he an ideal 
approach. Finite difference methods are ruled out by the 
infinite domain of the problem. Even if a maximum cell size 
can be fixed as for spherical cells and the integral term in 
(2.25) replaced by a finite approximating sum, the growth 
term has to be written using a formula from numerical differen- 
tiation. This may introduce large errors. 

The method of successive approximations would reciuire 
much more computational time and the final results would be 
in numerical form. It is difficult to predict whether the 
method of successive approximations would converge in the case 
of an infinite domain. 

Using MWR, a compact analytic form of the solution could 
be obtained with the constants in the solution related directly 
to the leading moments of the distribution function in question. 
If sufficiently accurate results are not available, it atleast 
may provide a reasonable initial approximation of the solution 
which could be improved by successive approximations. In 
either case, the saving in computation time would be signi- 
ficant. These factors indicate that MWR represents an attractive 
approach to the solution of the problem. 
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CHAPTER - IV 


SOLUTION OF POPULjTION BALANCE EQUATIONS 
LOR A MICROBIOL PROP^GiTOR USING MWR 


4.1 The Steady State Propagator - 

Por the steady state propagator, the mass distribution 
is given by equation (2. 20) 

CD 

[r(m) W(m)3 = 2 j p (m ! ) p(m, m 1 ) W(m') dm 1 

m 

- [p (m) + W(m) ...(2.20) 


The growth rate as shown earlier 
r(m) 


m 

e 


. .(2.23) 


Using (2.23), the fission probability function given by (2.1l) 


-( 


in - m, 


C .^2 


P' U, C s ) = 


r(m, C g ) 


G n/tt - / m “ m c 


erf c ( 


€ 


becomes for the steady state 


-(- 


m - m 


c_) 2 


n (m) 




1 

0 


m e 


erfc (g— £<0 


, , (4.1) 


The values of the cell parameters € , m c for rod-shaped cells 
are taken from the work of Eakman ( 7 ) end listed in .Appendix C. 
Por convenience in numerical computation, a change of 
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w(ex) 


When m f = m, y = x 
and as m — ?>co , x — 


..,(4*2) 


Equation (2.20) is accordingly written as 


§■ £xw(x)j = 2 1 P (G y) p(£x, £ y) w(y) dy 

J 

- jj 1 ( £ x) + 5 J w(x) 


fP(ex) 


. . (4-,. 3) 


Using (4.2) in Equations (4il) and (2,14)>(4*3) becomes 


5 fe [™ { U 


CD 

! ~ 2 r n / 


| j“i (y) P*(x, y) w(y) dy 


where 


^ (x) 


1 \fn 

and p ' (x, y) = j 


T(x) + £ 


e ~ (x - x c ) S 

erfc (x - x ') Q 


30 x 2 (v - x) 2 


Y = .— 

c e 


...( 4 . 4 ) 


..( 4 . 5 ) 


...(4.5a) 


Equation ( 4 . 4 ) is a linear integro-diff erential equation in 
the distribution w(x). Cells with zero mass cannot exist. 


Therefore, 


r(o) = 


... ( 4 . 6) 


Differentiating ( 4 . 4 ) once with respect to x and assuming that 
the second derivative of w(x) is bounded at x = o, it can be 
established that the derivative of w(x), 

w*(x) =0 at x = 0 ,..(4.7) 
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Note that p’(o, y) = o, as cells of zero mass cannot he 
formed at all. Leibnitz rule is made use of to differentiate 
the integral term in (4.4). It can he shown by a second 
differentiation of (4.4) that 

w"(x) 0 o at x = o ...( 4 . 8 ) 


.Assume that the distribution w(x) can he expanded in 
the form 


co 


w(x) * °n 0 n (x) 


. . .(4.9) 


n=o 


The domain of x is (0, co ) . The approximating functions 
0 n (x) are chosen as 


0 n (x) = e” x L^(x) ...(4*10} 

where L n (x) is the nth Laguerre polynomial. Nor a discussion 
on Laguerre polynomials, see -Appendix A. The family of 
functions ( 4 . 10 ) is complete in (0, cd ) ( 27) and an infinite 
linear combination as in (4.9) may be used to approximate 
w(x) in ( 0 , 00). In practice only a finite linear combination 
is possible. The trial solution therefore is 


N 

w*{x) - C n 0 n (x) 

n=o 


...(4.11) 


where N is the total number of functions 0 n used. 


Substitution of ( 4 . 11 ) into (4*4) gives 
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N 

0 


n 


L. 


n=o 


0 n '(x) § + 5 0 n (xQ 


CD 


N 


Z I p (y) p’(x, y) 


n=o 


c n v 


x 


IT 

21 C n 

n=o 


1 (x) + g-j 0 n (x) + R(x) 


R(x) is the residual of the integro differential equation (4*4) 

on using the approximation (4.I1) , 

letting 

0^(x) g + | 0 n (x) - 2 rp'(y) p’(x,y) 0 n (y)dy 

+ p (x) 0 n (x) x = Pn(x) ...( 4 . 12 ) 


R(x) can he written as 

R (x) = 21 °n 

n=o 


...( 4 . 13 ) 


The residual is averaged to zero using the MWR criterion, by 
weighting it with members of a complete set of functions AjK. , 
over the whole domain, 
co 

r 

j R(x) ,n ^-j c (x)dx =0 ...( 4 . 14 ) 

o 

for k = 1 , 2, N-2 
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Using (4.13), condition ( 4 . 14) is written 

<(R, 'f£(x)/> = Jr<^(x), Fn(xj> 0 n = 0 

n=o 

for k = 1 , 2, N-2 . . . (4.15) 


where 


<^u, v^> denotes ^ 

in future the inner product J u(x) v(x)dx 
Defining a^ n = ^^-(x) , Pn(x)^> , system (4.15) is written 


as 


N 

n=o 


kn °n = 0 


k = 1,2,, (N-2)...(4.16) 


giving (U-2) equations in (N+l) unknowns. The residual has been 

weighted with only (W-2) functions in order that the unknowns 

C n can he made to satisfy additional conditions already known 

regarding the distribution. It is expected that the rise of 

additional information about the solution will reduce the number 

of approximating functions required for convergence to the 

final solution. Condition ( 4 . 6 ) applied to ( 4 . 1 1 ) yields 

the (N-l)th equation in the C n 's 

N 

w*(o) = C n 0 n (o) = 0 ,.,(4.17) 

n=o 

Likewise condition (4.7) becomes 



n 


0 n * (0) = 0 


n=o 


...(4.18) 
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The (U +l)th eouation uses knowledge of the first moment of 
the distribution, C given by 2.1). 

The first moment of w(m) is 


co 


CO 


= C = / m V/(m)dm = e £ 


xW(x) dx 


...(4.19) 


Using (4.1l) 


C = 


co N 


2 r 


* J 


X 


G n ^n (x),Sx 


n=o 


2 R 
n=o 


°n 


...( 4 . 20 ) 


The orthonormal property of the approximating functions 
makes (4. 20) an exact condition regardless of the number of 
approximating functions used in (4.11). The first two Laguerre 
polynomials are 

L o = 1 

L.| = (l - x) 

The first moment of w(m) is 


A 


CD 

r 




w 


(m) dm 


CD 

& f w(x)dx 

J 

0 


...( 4 . 21 ) 


€ <w(&) , L q )> 

■§f °n <X* k> 


n= 0 
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Using 


X -x 
< e 


L , l\ 

m* ny 


= 0 
= 1 

^ o — p 

g“ ~ 0 


for m ^ n 
for m = n 


<X(x), = '—■ = C Q ...(4.22) 

<^w(x), L^(x)^> = ^w(x), (l-x)^> 

= <^w(x) , l)> - <(x, w(x)^> 

°n4: X 


X 

IT 


The left hand side 


n=o 


= C, 


The Right hand side 
(4.22). Therefore, 

°o " °1 = 


°o - 


C 


0 

e 2 


using the results (4.19) end 


**.(4.23) 


is the exact moment condition. 

Equations (4 . 16), (4.17), ( 4 . 18 ) and (4.23) form the 
system of algebraic equations to be solved for the C *s. It 
is interesting to note that the first constant C 0 multiplied 
by the constants, gives straightaway the first moment /i 
or the number of organisms in the propagator. The value of C Q 
as the number of trial functions is increased should approach 
a constant and will give an indication of the convergence of 
the approximation. 

The system to be solved is put in matrix form 
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<+i> v<^i ’ F i> 

f 

1 

I 

1 


" C o ~ 


0 

<+2’ ^<+2, ?!> 

<+> F a> 

<$^2 » F N> 

C 1 


0 

0 

<5 A N-2 ,i 2^ 

<^y_2 >p ^ > 

°N-3 

= 

0 

! * 

0 o (o) 0,,(o) 

0 2 (o) 

0jj(o) 

C U-2 


! * 

0^(0) 0’(o) 

02 * (o)~ - - 


°N-1 


• 

1 -1 

0 

0 

C N | 

f 

[ 

f X*r 

c 

£2 


Or 1 ~U = b ... (4*24) 


0 1 

h3 

il 

O 

O 

0 

— i. 

*0 

1 

! 

i 

i 

i 

• - , C ;J ) 


b = ( 0 , 0 , - - - - - 

■ - C/ e 2 ) 


a ±j = <$"i {x) > Vi (x ^ 

for i = 1 

to N -2 


3 = 1 

to N +1 (4.25) 

a N- 1 , j = 0 3 _ 1 (o) for 3 

= 1 to N +1 

...(4.26) 

a N, 3 = 0J_^°) for D 

= 1 to N +1 

. ..(4.27) 

a N+ 1 , 1 = 1 



a N+ 1 , 2 = _1 



a N+ 1 , 3 - 0 for j = 3 to 

(N+ 1 ) 

. ..(4.28) 


Equation (4.25) is expanded, using (4. 12 ) 

a ±3 = <6^i (x) - 0 j-i (x) §y + <i/'i (x) > 

y CO *V 

- 2^2 ^- i (x), f p'Cy) p'(x, y) 0 ; j_ 1 (y) dy ^ 
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The weighting functions to he used are of the form 

9^ j^(x) = e~^ x x^ - "* 0 ... (4.30) 

0 j_^(x) = e” x h^_ 1 (x) using (4-10) 

L^_^(x) is a (j-l) the order polynomial having real 

distinct roots in (o, 00). Each of the terms in (4.29), 
therefore corresponds to an integral having an oscillatory • 
function as the integrand. The evaluation of such integrals 
poses numerical difficulties in the use of ordinary numerical 
integration formulae such as the Newton-Cotes formula. This 
difficulty is overcome for the third term for example hy 
evaluating integrals of the form 
00 00 

j e“^ x x i “ 1 ^J y)e” y y m dy 

o x * . * (4*31 ) 

which has a well behaved function for the integrand, and adding 

up these terms as 

j 00 cd f 

JSE G d _ 1>k f e~ Xx x 1 -^ p'w e- y y ^dy 

k=1 0 x ...(4.32) 

where G- . , is the coefficient of x in the polynomial 

J“ * 

Vi (x) • 

Using (4.5a) for p'(x, y) and changing the order of integration, 
(4.32) becomes 
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00 


k=1 


G. - i 

0-1 f k 


P (y) 


-y k-1 

n v XT 


7 


J 


y dy 


y 

J e" ^ x x 1+1 (y-z) 2 dx 


O 


The inner integral can he evaluated analytically. The double 

integral was evaluated for each k by ? 32 point Gaussian 
quadrature formula (See -Appendix B) and checked to three 
decimal place accuracy with the evaluation of the same by a 
7 point Newton- Cotes integration formula. 

Similarly, the fourth term in (4-29) is expressed as a 
summation of integrals 
0 00 

2E G^_ 1 k ( e “^ +l)x p / ( x ) dx ...(4.33) 

k=1 ’ { 

A further advantage here of breaking up the function into 
several integrals is that only (2N-2) integrals corresponding 
to (2N-2) values of (i+k-l) need be evaluated as against 
(N-2) (N+l) integrals to be evaluated if the fourth term of 
( 4 . 29 ) is considered as a whole for each i and 3 . 


The first two terms of ( 4 . 29 ) for given by ( 4 . 30 ) 

can be evaluated exactly by additions of integrals of the form 
oo 



-(X+1) 


x 


m 


x 


m. ! 


c*+i) m+i 


...(4.34) 


o 

The matrix elements (4-15) are now known, Equation ( 4 . 24 ) 

was solved for N ranging from 4 to 18 using a matrix inversion 
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routine. Once the constants C n were known the approximate 
solution was generated by substitution into (4.11 )• The 
results are shown in Figures ( 4 - 1 ) and (4.2) and discussed in 
Section ( 4 * 6 ) . 


4.2 Steady State Propagator, Daughter Cells of Equal Massm /2 

v 

This model is described by Equation (2.3l), a functional 
differential equation 


_d 

dm 


jjr(m) W(m)J = 4p(2m) W(2m) - p (m) + 


W(m) 


...(2.31) 


The formulation of the problem for solution is exactly the 
same as for the integro-differential equation, except that the 
double- integral term in ( 4 . 29 ) for the matrix element a • • is 

-L J 

replaced by the integral 

4<^' 1 ( x b p'(2x) 0 ;j _ 1 (2x)/ ...(4.35) 

which is evaluated using Gaussian quadrature and checked by 
a Newton-Cotes 7 point rule. The conditions 
w(o) = 0 
w*( 0 ) =0 

and knowledge of the first moment <^x w(x)^/ > give the last 
three equations of the system ( 4 . 24 ) as before. The results 
obtained are plotted in Figure (4-3) and discussed in 
section ( 4 . 6 ). 
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T ABLE 4. 1 

VARIATION 0? ZEROTH MOMENT OF w(x) , 
C Q WITH' NUMBER OE TRIAL FUNCTIONS, 

N.WEIGHTIN& FUNCTIONS ARE OE THE 
FORM e - ^ x 11 


Number of 
trial 
functions 

N 

C 0 , Zeroth Moment x 10 ^ 

7\ = i 

A= 2 

A = 5 

4 

1.5940 

1 .5216 

1.5196 

5 

1.5225 

1.5509 

1.5085 

6 

1.4550 

1.4715 

1 .4848 

7 

1.4485 

1.4517 

1.4515 

8 

1,4160 

1.4441 

1 .4522 

9 

1.4158 

1.4026 

1.4555 

10 

1.4162 

1.5645 

1 .4902 

11 

1 .4150 

1.5900 

1.5585 

12 

1 .4152 

1 .4052 

1.5954 




TABLE 4.2 


COMPARISON OF SOLUTIONS TO THE STE4PY-STATE 
INTEGRO DIFPERENIIiL EQ UA TION (2.20) FOR NO . 
OFJPPROXIMiTING FUNCTION'S Nr 12, 


m 

cell mass 

_ ip 

gm x 10 

W(m) 

x 10-25 



Ratio 

Weighting Functions 

\= i X = 2 

e-t x x 11 

A_ = 3 

Pinal Solu- 
tion after 

2 successive 
approx, of XfI 

L.H.S. 

R.H.S. 

f or A = 2 

0 

0 

0 

0 

0 

1.0000 

0.4 

0.5449 

0.4612 

0.4512 

0.4854 

1.0002 

0.8 

1.3456 

1.4068 

1 -4075 

1.4861 

0.9995 

1.2 

2.2905 

2.3045 

2.2945 

2.4496 

0.9974 

1 .6 

3.1605 

2.8687 

2.9121 

3.0303 

0.9210 

2.0 

3.1845 

2.9277 

3.0572 

3.0955 

1.0328 

2.4 

2.4149 

2.4391 

2.5444 

2.6578 

0.6393 

72.8 

1.4142 

1.6213 

1*6067 

1.5951 

1.1511 

3.2 

0.6205 

0*8085 

0.6987 

0.3829 

1.7632 

3.6 

0.1639 

0. 2285 

0.1152 

0.1494 

1.7631 
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4.3 Veri f icat io n of Results Ob tained by MY/R by Successive • 
■Approximations ° 

a) Integro- Differential Equation : Equation (4*4) is 
written as 

co 

|f + | pf +f (x) J w(x) = y p (y) P ! (x,y) w(y)dy 

x ...(4.36) 

N 

I3ae solution w°(x) = C 0_(x) ...(4.37) 

-fiC II XI 


11=0 


obtained from MWR in Section (4.2) is considered a zeroth 
approximation to the true solution. The first approximation 
wPx) is obtained by substituting (4.37) in the right hand 
side of (4.36) and solving 


d\P , Q 

dx + X 


jj§ + p (*0 wpx) = ~ J p 7 (y) p * (x,y) w°(y) dy 

x ...(4.38) 

Equation (4.38) is an ordinary differential equation to which 
the solution is given by 

f" Q rfe ! ) dx * 

n 3C " Tr 


(x) = 


x 

r 


J 


2Q x n e J 


dx" J p / (y)p(x",y) w°(y)dy 

...(4.39) 


Hie integration constant equals zero since w(o) _=.o 
Eor p / (x) given by (4.5) 


x 

r 


Q r (x 1 ) dx~ 


x' 


Perfc (x - x ) 

ln effcTaP^TT 


...(4.40) 


and equation (4.39) for the (n+l)th approximation, given the 
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nth approximation, using p’(x, y) from (4.5a) becomes 


x ■, cp 

^(x) erfc (x - X 0 ) f ^ j Jp'iy) (y-x") V(y)dy 

vj 7\ X J C 


C ?! 
x” 


. . .( 4 . 41 ) 


The integration of the outer integral in ( 4 . 41 ) is done using 
an integration subroutine which on a single integration gives 
values of w n+ ”*(x) for different values of x. 


Two successive approximations were found to be enough 
for convsrgence to the final solution within 3 decimal place 
accuracy. 


b) Functional Differential Equation * 

The solution obtained by MWR for the functional differential 
equation in section (4.3) was improved upon by 


w n+1 (x) = %~z . erfc( 

x tJ/x 


f -5^-4 w n (2x") ax' 


x z c 1 j erfc (.x" -x c ) 


0 ...( 4 . 42 ) 

derived similar to (4.41 ) for the integro differential 
equation. The zeroth approximation assumed the value zero 
for W(m) for the region 0 to 1 x 10 where negative values 
of W(m) were generated by MWR, 


4.4 of Approximating Functions 0 n (x) = e -x x n ”. 

The approximating functions originally selected were 
0 n = e~ x x 11 , which are the basic functions reauired to be 

orthogonal ized in order to generate the Laguerre polynomials. 

# 
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. "“X 21. 

With 0 n = e x the first moment c ondition is 

c° „ 

5E C n 0 n (x £> = A 

n=o ^ 

For a finite number of approximating functions the above 
condition is satisfied only approximately, becoming more and 
more exact as the number of trial functions is increased. So 
in effect, the use of approximating functions which are not 
orthogonal to each other results in both the differential 
equation and the first moment equation being satisfied only 
in an approximate sense. In such a Case, thore approximating 
functions are needed to generate a reasonably accurate solution. 

The results obtained bear this out. There was no definite 
trend in the improvement of the solution as the number of trial 
functions was increased. With the use of orthogonal polynomials 
as approximating functions, the moment condition became an 
exact one and the solution showed better convergence properties. 
This in addition to the fact that the moments of the distri- 
bution are directly related to the constants of the trial solution 
formed the basis of selection of orthogonal polynomials as . 

W 1 

approximating functions. 

4 .5 Use of Collocation Method for Integral Equation s 

A collocation method' in which the residual (4*13) was 
equated to zero at N points in the domain was tried first 
with 0 (x) = e -x x n as approximating functions and then with 
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0 n (x) = e~ x L n (x) . The points were chosen 

a.) at regular intervals 

h) as zeroes of the Nth Daguerre polynomial when 
using 0 n (x) = e“ x D n (x) , n = o to N as suggested hy Stewart ( 26 ) o 

The method failed to work in "both the cases, giving 
erroneous results even when 25 approximating functions were 
used. It appears that the Collocation method is not applicable 
for this problem with a reasonable number of trial functions 
and a more elaborate procedure of averaging the error over the 
whole domain is required. 

4«6 Discussion of Results for Steady State Solutions * 
a) Integro- Differential Equation ? 

Three different weighting functions were considered 
first, of the type ( 4 » 30 ) with A. =1 , 2 and 3 . Figure ( 4 - 1 ) 
shows the variation of W(m) for different number of trial 
functions N for X= 2 and 6 - 2 hrs» It was observed that the 
solution settled down to its final values, differing only in 
the third decimal place with about eleven or twelve approxi- 
mating functions. Table ( 4 . 1 ) shows the variation of the 
zeroth moment of W(m) given by eC o , with different N for 
< X.= 1> 2 and 3. Figure (4.2) and Table (4.2) compare the 
solutions obtained for ,V= 1, 2 and 3 with the number of trial 
functions fixed. From Table (4 * 2 ) it is seen that the solutions 
using weighting function e x and e x are closer to 
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N = 12 FNS. 

0 (X)=€ x Un(X) 

WEIGHTING FUNCTIONS 

4 . e x Ln(X) 

2. € X Xn 


3. e 2x Xn 


m, CELL MASS, ($rn) x ICT* 

FIG. 4-2 SOLUTIONS FOR STEADY STATE PROPAGATOR 
COMPARISON OF WEIGHTING FUNCTIONS 




the final solution than that using e” x x n as a weighting 
function, especially in tie rising portion of the solution. 

Figure (4*2) includes the curve obtained by using 
Galerkin's method where the weighting function ^ * is the 
same as the approximating function 0. . In this case, both 

J 

are polynomials and the integrals to be evaluated involve 
integrands with i+j roots, for a product . 0j). Thus 

a double summation of integrals of the form (4-31 ) and (4.33) 
is required and this was more in error as even a single 
summation process like (4.32) was found to involve addition 
and subtraction of nearly eciuel quantities. So one expects 
more deviation from the true solution for this case. These 
difficulties indicate that Galerkin’s method is less accurate 
where oscillating polynomials are used as approximating 
functions. The error involved in the summation process increase 
as more trial functions are used because of the higher order 
polynomials so introduced. In such a situation increasing the 
number of trial functions may not always give an improved 
solution as it was observed with the Galerkin's method. In 
this respect, the method of moments seems superior to the 
Galerkin's method. 

Besides observing how the solution varies with the 
number of trial functions and different weighting functions, 
an index of the extent to which the differential equation is 
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substituting (4.11) in (4.4) and talcing the ratio of left 
hand side to right hand side. Table ( 4 . 2 ) lists these ratios 
at different points in the region. For values of mass upto 
about 2.4 x 10” gm (9 little beyond the peak) the ratio is 
almost one, showing that the solution is very close to the 
true curve. Beyond this the ratio varies widely. 

The final solution is obtained in as few as two approxi- 
mations. The zeroth approximation from MWR deviates at the 
most by 5 % from the final solution. This shows that the 
approximate solution is itself valuable considering that ex- 
perimental results are bound to be in error atleast by this 
amount . 

The index, ratio of the left hand side to the right 

hand side is not a good indicator of the correctness of the 

— 1 2 

solution, especially in the region of mass beyond 2.4 x 10 gm. 
When substituting the approximate solution ( 4 . 11 ) an (4.3) 
and taking the ratio of left hand to right hand side at a 
particular value x^ of x, the integral term in ( 4 - 3 ) considers 
values of w(x) for x from x^ to co . The percentage error 
involved in evaluating the integral increases as x^ becomes 
larger as w(x) for large x^ is more in error than for small 
Xj_. Table ( 4 . 2 ) reveals that the difference between the final 
solution and that of MWR for cell masses beyond 2.4 x 10 is 
not appreciable enough ( 5 % deviation) to account for the 
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Figure (4.5) shows the variation of the distribution 
¥(m) with holding time 9. Decreasing the holding time 
increases the number density of cells throughout the mass 
range. This is because the growth rate r(m) = m /Q and the 
fission probability at steady state (4.1) increase with 
decrease in 9, resulting in more cells of all masses. This 
increase more than offsets the loss due to the increase in 
the washout rate term I/O. 

b) Functional Differential Equation ’ 

The situation with respect to the functional differential 

equation is some what different. The solution generated by 

MWR is in great error in the region in the region of mass 

0 to 1 x ICf* 12 gm. However, six successive approximations are 

sufficient to make the solution settle down to its final value. 

The final solution (Fig. 4 . 3) shows that W(m) is very small 

(as low as 10 organisms per litre per gram) for values of 

mass from 0 to 1 x 10~ 12 gm. The maximum value of W(m) is 

-1 2 

108 times greater than the value at m = 1 1 10 gm as 
contrasted with the solution of the integral equation where 
it is only four times the value. The nature of the solution 
compares well with that obtained by Eakma,n(8). 

One may conclude from this that it is difficult for MWR 
to generate such a, solution peculiar to the functional 
differential equation, it the most it may provide a rough 

n f the solution in some regions. it may be difficult 


eo 



12 

m, CELL MASS Cgm) x 10 


FIG. 43 METHOD OF SUCCESSIVE APPROXIMATIONS FOR THE 
FUNCTIONAL DIFFERENTIAL EQUATION 



Method of successive approximations for steady state integr 

differential equation. 


2 


3 


4 


(m),CELL MASS GM xio ' 2 

STEADY STATE CELL MASS DISTRIBUTION FOR DIFFERENT 
HOLDING TIMES 
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for a linepr combination of a smell number of polynomials like 
(4*1l) "to generate p solution of such wide variations in its 
mpgnitude . 

c) Compprison of Models : 

Figure ( 4 . 6 ) compares the results for the three models 
described in section(2*5) by Equations (2.28), <2.31 ) end 
(2.2l). The general model IV ( 2 . 2 1 ) gives results close to 
those of model III in the region of smell mess. This is to 
be expected since P (m) is very smell for smell m Pig. ( 4 . 7). 

| * (m) is zero for model HI till m = m Q where it is infinite . 
It does have a finite value end this accounts for the curve 
for model IV being Ho ve thet of model III in the regions of 
smell mess. 

* 

Pig. ( 4 . 7 ) shows thet p (m) storts rising steeply after 
-12 

m = 2 x 10 gms. This indicates that very few of the dividing 

-12 

cells have masses below 2 x 10 gms. Thus if one assumes 
that cells divide exactly into two h-elves, there will be very 
few cells with masses below 1 x 10 gms as shown by the 
final solution to the function differential equation (2,31 ) 

The integro-differential equation considers a distri- 
bution of daughter cell masses and so one can expect more cells 
of small mass in this model. 



M. Eqn.(2-21) 



00 

CM 


CM 



oiinamisia * (^)m 


CELL MASS (gm ) X 10 
Fiq.4-6 
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4.7 The Unsteady State Propagator . 

The partial integro-differential equation describing 
the mass distribution of cells in the unsteady state (2.20) is 


+ ^ L r(m ’ °s ) w(m « t) _ 


Co 


r(m, C g ) 


'2 K c s (t) 


= 2jp(m',C g ) 

m 

_fp(m, C B ) +l^W(m, t) 

L ,..(2. 20 ' 


h k s +-sy t ) 

f ( t) m 


A c m 


...(2.15) 


Changing the independent variable according to 
m = 6 x, m f = £ y 
(2.20) is written 

CD 

%% + f(t) — fx w(x, tQ= 2 f p ty, t) p 1 (x,y) w(y,t) dy 

O v o x * — J 


X 


-fp'U, t) + 5 

¥(m, t) = W(£x, t) = w(x, t) 
p’(x, y) is given by (4.5a) 

2 


w 


(x» t) ...(4*45. 


2 e- (x " x c } 


V* <c \ffr erfc be - x c 7 


f(t) 


X 


. ..(4.44) 


ra 


On using (2.15) and - x_ in 


in (2*1l) 


\SL 

bt 


c 

£ * c 

Equation (4.43) is re-arranged in the form 

co 


f(t) 


x 






4" 


2 J^p r (y» t) p ’ (x,y) w(y , t) dy 
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where L is an integro-differential operator such that 


v 00 


Lv = - f(t) x + 2 

o -*■ 


p Cy, t) p'(x,y) v(y) dy 
x 



- £pfx, t) + £ +-rjf(t)3 v(x) 

...(4.46) 

Assume that the 

solution to (4.43) is of the form 

TVT 


W*(Xj, t) = g ( x) 

N 

+ S. °n (t) 

.. .(4.47) 


n=o 


where w* (x, o) 

= g(x) 

.. .(4.48) 


with C n (o) = 0 for n = o to I 

0 n are the orthonormal system of Laguerre polynomials as before. 


Ihe selection of the form of (4.47) merits attention. 

in expression of the type ]>[ simplifies the selection 

of trial functions. Solutions based on 

w*(x, t) = c n (t) 0 n (x) 
n=o 


can resolve the time dependent behaviour more completely since 
the errors can be orthogonalized. at every instant of time insteai 
of being orthogonalized over a finite period of time. ALso in 
most problems a smaller set of functions may be used because 
the 0 n 's do not depend on t. 

Condition (4.48) ensures that the initial residual 

11 

w*(x, o) - g(x) - C^( ° ) 0 n (x) 

n»o 


is exact"! v 7;er0. 
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Substituting the tripl solution w*(x, t) in (4.45) 


± c- ( t) « a (x) = l r S (x) + si 

n=o . n=o 


C n (t) 0 n (x) 


The residual R(x, t) of the equation (4.45) is 


C*(t) 0 n (x) - Ig(x) 


0 n (t) L0 n (x) 


.(4.4-9) 


The MWR criterion 


OO 

r 

J R(x, t) '(/'% ( x) dx » 0 


...(4.50) 


j = 1 to H-1 


of weighting the residual with members of a complete set of 


functions ^ , gives the sys 


xn V 

<0 n /j> = °n 

X ^ n=0 ...(4.51) 


j = 1 to N-1 

Weighting has been done with (R-l) functions only, so that the 


known condition 


(4.52) 


w( 0 , t) = o for all t 

and knowledge of the first moment of £he distribution from 


(2. 18) and (2.19) may be used. 

rt is to be noted that a condition corresponding to ( 4 . 7 ) 

regarding the derivative of the distribution respect to x at 
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x — o does not exist in the unsteady state. Differentiating 
(4«43) partially with respect to x and assuming that the second 
partial derivative at x = 0 is bounded, the following 


condition is arrived at 
w ! 




t 2 f(t > + 5 J - 


~Vw 


~b t 


x=o 


...(4.53) 


x=o 


The right hand side of (4.53) need not in general be zero, 
"c> w 




is therefore not zero at x = 0 . 


Condition (4.52) provides the Nth equation in C. 


n 


g(o) + 


N 


0 n (t) 0 n (o) = 0 


n=o 


or in terms of derivatives with resnect to t of C ’ s 


N 


C'(t) 0_(o) = 0 


« . * (4 . 5-r ' 


n=o 


The (N+l)th equation is derived making use of the moment 


condition 


co 

r 


w 


(x, t) .1. dx = 


A> (t) 


• . * (4-f 55: 


w* 

0 


/^(t) is the number density of organisms in the propagator 
at time t. Using (4.47) for w(x,t) in the left hand side of 


(4.55) 


CD 


} w(x, 


= <(g(x)j}> + 21 c n ^n L o/> 


t) L Q (x)dx 
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by orthogonality property of L Similarly, 


J w(x,t) i 1 (x)dx = <w(x, t) p (l - x )^> 

° _ /^o(t) _ 

6 p 2 


... 4 . 


ilso 


uu 

w(x, t) (l - x) dx = <^g(x) P (l- x)^> 
o N s 

+ SI C n<A^ 


Thus 


>k(tl J-^|I . </( x) , (1 _ x )> + Cl (t) 


Using (4-56) for 


■X 0 

~v— equation (.4.58; gives 


0 Q ( t) - C.j(t) = - <^x, g(x^> + 


...(4.5 ; 


...(4.: 


C(t) is the first moment/^ ( t) and represents the biomass con- 


centration. 


Differentiating (4.59) with resnect to t 

eye - c^tt) = i, 

is the (N+l)th equation to be used. 


. ..(4.60 


is obtained at any time from the simultaneous solution 
of (2.18) and (2.19) 


iL + (2 

Q '‘ITp 
<°s - °s> 


0 C 

s 


s s 


/V 0 


ffL 


0 C C 
s 


s s 


...( 2 . 1 - 


...(?, 
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Equations U.5l), (4.54), ( 4 . 6 O), (2.18) and (2.19) are to 
be solved simultaneously for C(t), C (t) and the F + 1 
variables C Qf (t) , . . . C^( t) . The forms of the terms in the 

system of equations (4.51 ) will be examined. The second term 
on the right hand side is. got by summing integrals of the type 

<6 0 n (x) / V' 3 > 

= - f(t) <4 X 0 n <( x ),'/'.(x) y 

2 J n'(j, t) p'(x, y) 0 n (y) dy^> 

-<^p(x, t) 0 n (x) , 'VK.(x)^> 

" I <^ 0 n (x) ’ ^ (x) ^ 

- f(t) <Q^(x) , (x)*^> ...(4.61) 

The first term on the right hand side is 
L g (x) 

which is of the same form as ( 4 . 6 1 ) with g(x) replacing 0 n (x) . 
Examination of ( 4 . '6 1 ) reveals that the fong-^f the integrals 
to be evaluated numerically corresponds to the type given by 

-v 

(4.32) if one uses weighting functions of the form e” x n . 
This is so because p*(y, t) is expressible as (4.44) where 
the function of time can be separated out and the remaining 
portion is a function of x only. Thus there is no need to 
evaluate integrals other than that used for the steady state 


solution. 
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Define. 


•(x - x) 2 


M(x) = g- 7 7 v 

erfc lx - xj 


2 


« « • (4*62) 


axid 


co 


H n (x) = x 0 n '(x) - 2 J M(y)yp ! (x,y) 0 n (y)dy 

x 

+ x$(x) 0 n (x) + 0 n (x) 

The second term on the R.H.S. of (4.5l) is now 
N 


Z C n|_< H nW-tj(x>} - jjrJT C n <^0 n ,t>? 

3 = 1 to 1-1 ... ( 4 . 63 ) 


n=o 


The system of equations (4- 51 ) (4.54) and ( 4 . 6 O) is express" 
as a matrix differential equation, using (4.63). 


I £5 

A ^ 


f(t) BC' ^ DC + f(t) 


w 


Q 


v + t 


...U-£4) 


where C = (C Q (t) , C^(t)- - - - - - C^(t)) 



<?2 ^ I s ) - * ' i)> 

<^2> <*1 & z 7 



^N-2> 

<^2;^ S=6/ > ’ N-^> 

-Si 

0 

0 

^Sl 

0 

0 2 (0) 0 N (o) 

1 -1 

0 0 
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D is a matrix with elements 


d, » • " — * ■ 

10 10 


=0 

10 


for i = 1 to H-1 

0 = 1 to N+1 

for i = IT, N + 1 

1 = 1 to IT + 1 


^(x).^^ 4-' 4 n U),t^l 

<H 0 (x),^ 2 (x^> <H 1 (x),f' 2 (x)> <H n (x),^(xj>| 


4o^,T%. 2 ^)> ^1 (x) <^n ( ^fk.2^)> 


!The general -term for i = 1 to N-2,j = 1 to N+1 , ( 

requires only the integrals used in the steady state problem, 
w is a vector such that 
w fc = - /g'(x)x, 0\(x)^> 

^ m 


.(x) • 2 


yM(y)p’(x,y) g(y)dy 


X 

<^x M(x) g(x)^^ k (x)y> - ^g(x)^ k (x)^> ..'.(4.65) 


for k = 1 to N- 1 . 

W 1T = W N+1 = 0 
V is a. vector such that 

= S g(x) 'V'V(x)')> for k = 1 to IT - 1 


. . . ( 4 . 66 ) 
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t is a vector such that 

\ = 0 for k = 1 to N 

+ I_ dC(t) 

N+1 ~ £ 2 &T~ 

G(t) is the biomass concentration at time t. 

Multiplying (4.64) throughout by 1 ~ 1 

= - f(t) I ~ 1 1 c - ^ 1 - 1 5 3 
+ f(t) I -1 w - A I' 1 v + 1 - 1 1 ...(4.67) 

which arc a set of first order ordinary differential 

equations to bo solved for C simultaneously with- (2.18) and (2*10 

which 'determine C(t), C (t) and hence f(t) in (4.67). 

s c. 

The initial conditions are C(o) = 0 and ' - 

co 

C(o) = J g(m) dm = 6 <^g(x) J) 1^> 

0 

and C (o) is specified arbitrarily. 

The system of equations (4.67) together with (2.18) 
and (2.19) was integrated using a Runge-Kuttp.-GJ.il method. 

The solutions were obtained for the following types of initial 
distributions^ 

a) g(x) = k x e"~ X . 10 2 ^, a population consisting of pre- 
dominantly cells of small mass. Two values of k were used. 

The results arc presented in figures ( 4 .IO), (4.1l)> (4.13) 
and ( 4 , - 14 ) for approximating functions 0 n = e x I< n (x) and 

J J_ • « -C* 4-V ^ 4 *TTT-,rt ( Trh Q TT^- 


75 


b) g(x) k x e , 10 a population consisting of 
predominantly cells of large mass. Hie solution is shown in 
Figure (4-12) 

Figure (43 ) compares the solutions obtained at a 
particular time as the number of trial functions is increased. 
Figure (4-9 ) shows the effect of using different weighting 

functions, for the same number of approximating functions. 

Hie results are discussed in section (4.8). 


4-8 The Batch Propagator : 

a) Changing Environment : The cell mass distribution of 
batch culture in a changing environment was obtained by 
solution of equation (2. 32 ) by MWR. The formulation of t v 
method of solution is the same as for the continuous propagator 
equations but for the omission of the ( 1 /©) term in (2.20; 

(2. 18). The system of differential equations to be solved row 
is 


i dC 
A Ft 


f (t) "B D + f(t)w+t 


.. .(4.68) 


f(t) T 9 K °s (t) -A 
f(t) -L 2ir p 


as before 


dC f 0 ^m 

3t -p xy 



+TT) 



( 2-33 ) 




C s C 

xr+o 

s s' 


(2.54 ) 
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subject to 



C(o) = J g(m)dm where 

o 

w(m, o) = g(m) an initial distribution on and C (o) is 
an arbitrary initial substrate concentration. 

The results obtained are presented in Figure (4. 16). 

b) Exponential Batch Growth : 

Considering the environment to be essentially constant 
the distribution of cell masses is obtained as in (a ) , but 
omitting the substrate balance (2.33) and assuming C Q a 
constant in (2.34). The exponential growth phase is shown 
by the straight line portion of the plot of 'In N against time 
in Figure (4. 18). Exponential batch growth mass distribution 
can be written as 

W(t, m) = E 0 exp ( P t) g(m) ...( 2.35) 

Knowledge of fro.nn the slope of the straight line portion of 
Figure (4. 18) allows calculation of g(m) from w(m, t) . 

■K - 

Equation (2.36) for g(m) was solved independently by 

M¥R similar to the solution of the steady state propagator 

00 , 

balance (2.21 ) but using f g*(m)dm = 1 instead of the momen • 

condition (4.25) and the value of obtained earlier. The 
solutions for g4m) by the two methods were in good agreement 
(Fig. 4.19). initial distribution of the type W(x, o)= xe 

was assumed. 
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4*9 Discussion of Results for the Unsteady State Equation ; 

a.) Test of Corr e ctness of the Soluti ons * 

Whatever the type of initial distribution chosen, the 
unsteady state solution approached the steady state solution 
determined independently, to a good degree of accuracy. This 
is only a necessary but not sufficient condition for the 
correctness of the solution. However, several factors, 
discussed below, provide reasonable bases for assuming that 
the solutions were correct. These factors relate to, 

1. Doting the rate of convergence of the solution at 
any fixed time co-ordinate by increasing the number of apprext 
mating functions. .About ten functions were found enough to 
give a. solution which did not change appreciably on further 
addition of trial functions. This conclusion applies alsa 

to convergence of the first moment of W(m, t) to its final 
value . 

2. Solving the problem for different weighting functions 
Three different weighting functions were used. Dor the same 
number of approximating func bions, the different weighting 
functions gave reasonably close solutions at any time* The 
variation between the results was at most 

3. Solving the problem for a situation where the 
of the solution is known definitely. Batch growth under 
constant environment was treated as a special case of the 
general problem. The solutions obtained fitted very well with 
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the exponential phase mass distribution W(m, t) = N 0 e g*tm5 . 
On using the value of obtained from this fit, in the integro- 
differential equation for g*(m) and solving it, the distri- 
bution obtained corresponded to that got by the relation 

g*(m) = 11 „e~ j t 

0 

b) Interpretation of the Results for a Continuous Propagator 
1. Figure (4.1 1) shows that for an initial distribution 
having comparatively larger number of cells of small mass, the 
total number density in the propagator N(t) decreases first. 

This can be interpreted by saying that the probability of small 
cells dividing is small and thus in the initial stages the loss 
in number by washout is significant compared to the term account 
ing for gain by division. Starting with zero substrate in the 
reactor, one finds that growth is also slow initially and the 
substrate concentration a,t the outlet of the reactor rises to 
a maximum value. By this time the cells have grown enough to 
start multiplying by division. Rapid increase in the number of 
cells Increase the rate of consumption of substrate which startj 
dropping down rapidly. 

The distribution (Fig. 4.10) follows the behaviour of th< 
first moment which is determined independently here, il though 
initially the number density decreases, the first moment does 
of the formation of more cells of large mass by 


not, because 



79 


the propagator as well as the biomass concentration is brough 4- 
about by division of larger cells. As a result, substrate is 
consumed very rapidly and the biomass concentration and number 
density overshoot the steady state value. The steady state 
is reached in about three hours. 

2. in initial cell mass distribution of very few cells 
mostly of small mass takes a much longer time to assume the 
steady state distribution. (Figures (4. 13) and (4.14)). The 
number of density and mass distribution change with time 
similar to case (l) but the initial drop in the number of cells 
is not as significant as in (l) because there are very few 
cells totally. 

3. in initial distribution with cells of relatively 
large mass changes in time as shown in Figures (4.12) and u 
The number of cells of smaller mass starts increasing as a res'; 
of division of large cell and the entire curve shifts towards 
the left to the steady state value. .Accordingly the total 
number of cells F rises steeply at first and settles down 

to its final value. The steady state is reached much faster 
compared to cases (l) and (2) since the initial distribution 
is closer to the steady state distribution than case (l) or 
case (2). 
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c) Results for the Batch Propagator : 

1 . Changing Environment : 

Ihe mass distribution for a batch propagator is shown 
in figure ( 4 . 17 ). The biomass concentration c, the substrate 
concentration C g and the number density N, vary with time 
as shown in figure ( 4.1 0 . It is seen that after about 0.5 
hrs the substrate concentration decreases to very small values, 
the point at which the biomass concentration starts dropping 
rapidly. The growth rate becomes negative in the particular 
model used. The transition probability of cell division as 
in (2.1l) would then become negative which is impossible 
physically. It can be shown (32) that 

r (m, C g ) 

should be replaced by jr(m, C g )jin ^he expression (2.1l). 

A negative growth rate would reduce the number of 
cells of large mass, and increase the number of cells of 
small mass (see figure ( 4 . 16)). The total number of cells 
eventually becomes constant as no death has been assumed and It 
is more difficult for smaller cells to divide. 

2. Constant Environment : 

The solutions for constant environment fit very well 
with the expected form of number increase R(t) — e . 

After an initial adjustment period the plot of log vs time 
(fig. 4.18) is a straightline representing exponential growth, 



81 


in exponential growth. The asymptotic mass distribution g(m) 
calculated from W( m , t) on solution of (2.32) directly with 
C s constant is nearly identical to that obtained by solving 
the integro-differential equation (2.36) for g(m) with 
9 = 1.282 (Figure (4.19)). 

d) Numerical Difficulties Involved in the Solution of 
System (4.67) ^ 

For the approximating function family 0 n (x) = e~ x I (x- 
and weighting function family (x) = e"" x x n it was found 
that the system of equations became difficult to solve for 
matrices larger than order 7, The matrix 1 was such that its 
inverse contained large elements and extremely fine interval 
were required in the Runge-Kutta method for the solution to 
remain stable. This trouble did not arise for the systems 
using 0 n (x) = e~ x 1 (x) and weighting functions / 'f^ n (x)= a 
and e~^ x x 11 upto a matrix order of 11. Fortunately it was 
found that eleven functions were enough for convergence to tu 
final solution. The higher order matrices were found to have 
very small values for their determinants and indicate ill- 
conditioned matrices. Is a result, even after double pree: 
computation, the inverse was likely in error and probably 
caused instability of the differential equation system (4.6, 


for large step sizes. 
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e ) Computation Time ; 

The total time used for the problem was about 20 hou^ 
on the IBM 7044 computer. Most of the time was taken up in 
attempting to calculate accurately the double and single 
integrals involved in the matrix elements of the final systc _ 
to be solved. Double integrals evaluated by Gaussian quadra tu 
(See .Appendix B) take about 5-10 seconds each to compute as 
compared to 20-30 seconds each when using a Newton-Cotes 
formula. The single integrals take about half the time reqv" 
for double integrals. 

The integration of the differential eouations for a 
matrix order of 10 required about 5 minutes for getting vali''- 
of the distribution in step sizes of 0.005 hrs. from 0 hr. to 
5 hrs. 

The time required to solve the system of algebraic 
equations in the steady state case is not significant. 

4.10 Suggestions for Future Work 0 . 

The solutions to the steady state equations were checked 
quite easily by the method of successive approximations. The 
method of successive approximation could be applied for the 
unsteady state equation using the zeroth approximation from 
MWR. A partial differential equation would have to be solve' 5 
at each stage by the method of characteristics. 
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The problem could be solved for the case of spherical 
cells also which are commonly encountered. 

The method of weighted residuals could be extended to 
solve the system of vector differential equations describing 
the population balances for a structured, distributed model 
of a biological culture which is more general in its approach 

Partial integrb-diff erential equations of the type 
solved in this work arise in many other situations involving 
number balances. MWR could be tried on non-linear equation^ 
such as those encountered in a dispersed phase system where 
coalescence phenomena care considered. 
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LAGUJSfiRB POLYNOMIALS 


M important class of orthogonal polynomials which are 
complete in the space (o, co ) are the Laguerre polynomials, (27, 
28, 31), defined by the formula 

L n (x,e C, ft ) = e° < ” x x 1 “ v " df; / - c£x n+^-1\ 

n x } 


oC > O 


dr' 
f> > o 


n = o, 1, 2 . . (A. l) 


Por cC - — 1 one has the ordinary Laguerre polynomials. 


,n 


V x) = « , n 

dx 


x d ‘ (e~ x x n ) 


n 


> 


>=o 

For the genernl cose 


( 5 } fr ( - x) 


9 


...{A. 2) 


n x > 

L n (x ; cjc_, A ) = x ( n ? ) (p+9; 1; n -j)) (-oCx) ...(-fi.3) 

J>=0 ^ 


where 

(m:d: } ) = m(m+d) (m +2d) m +(<P-l)d 

>) = 1, 2, 

(m:d:o) = 1 

The first few laguerre polynomials are thus 


L 0 (x,aC , £ ) - 1 

(x,ct,^S ) = p i-aCx 

(x,oc,6 ) = £(/5+l) - 2(y6+l)<x +.^ x 2 


U.4) 


* • » 
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The Lngutrro polynomial \U, X, p ) has n distinct real roots 
in (o, co ) and any function multiplied by L n (x,oC,£) will show 
oscillatory behaviour. 

The family (i.3) satisfies the orthogonality condition 


oo 


j x‘ M i n (x,oc ,p ) i m (x,ot,|& ) ax 


0 


m 


ja J ,P ( /6 + n) 


cC 




^ n 

m ~ n • . . { A. 5 } 


To normalise the polynomials (A. 3) they should be divided by 


nl p(/3 +n) 

‘ ** 

Lngucrre polynomials satisfy the recurrence relation 
"^n+2 - 2n - f3 ~ 2) T n+ -| (x, oO ,^2> ) 

+ (n+l) (n+ ; 5) h n (x,0C r ^3) = 0 

n=o,1,2 .»» * . . ( A*6) 

and tho Lngucrre differential equation 

x L n " (x,oG ,f> ) + (p- exx) L n T (x,oe,/i ) 

+ n c?C L n (x, oC, jb ) =0 
n = 0, 1, 2 ••• . . . ( A. 7) 


It may bo shown using the orthogonality property that 


co 

J o-** x JB+ ^- 1 V x > oC ’( S } 
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Ex pan si on of Functions in Series o f haguerre Polynomials (28) 

Uuc of the most important properties of the haguerre 
polynomials, in the foot that a reel function f(x) defined in 
the infinite interval (o, <n ) con be expended in a series of the 
form 

r (x) = ZEL C n L n 
n=o 

0 C x < fX> ...(1.9) 

provided f(x) satisfies cv rtsin conditions. The coefficients 
C con h>; determined using the orthogonality property (1.5). 

The extension (1.9) is valid if f(x) is piebewise smooth in 
every finite intorvnl , x 9 J and suitably well behaved near 
the points x = o and x = co . 
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c/pggian QUADRA tube usin g daguerrei polynomials (50) 


i''or tht evaluation of integrals over infinite intervals, 
then are two approaches (l) Uee a knowledge of the integrand 
to bound the magnitude of the integral from some finite value 
to infinity by a positive constant 6>0 and then use a quadra- 
tun formal'" for the remaining finite- interval. (2) Use a 
nundrat.uro formula especially developed for the infinite 
interval. The latter approach is more convenient and is eco- 
nomical in terms of computer time. The former can be used as a 
check when necessary. 

The L -"guerre- Gauss quadrature formula is given by (30) 
co 

- . n 

j e~ x f(x)dx = yr H _ f ( a _) + E . . . (B. l) 

3 3 

O i 


WhL IT 1 t-lii, r • T s *rre the zorors of Laguerre polynomial 

of nth order. (See Appendix A) the H_.’s being given 


J H 


(nih 

u 1 ^ vi Ci 


.. .(b.2) 


The error term E is 

a (n!) 2 jl fzia L, (B. 3 ) 

" - C 2nTi 

^ lies in (0 , oo ) 

If f(x) is continuous, any desired degree of accuracy is 
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obi?- in- 1 .!-' by using a sufficiently high order formula, by 
ine rt u, ibo number of Quadrature points in (B.l). 

Using serjiH-nco o-f Gaussian formulae like (B.l) with increasin 
n will la 'ui io a canvergent seouence of approximations. 

Tb«, a.V and in’s in (B.l) have boon listed for n=2 to 32 in 

t) u 

tlu api'f.'idix of (29). 
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APPENDIX - C 


D ATA POR ROD- SHAPED CELLS (j) 

— 12 

m_, mean division mass = 3 x 10 gm 
t£ , measure of spread in division mass distribution = 

3 j 2 x lO' 1 ^ gm. 

0 , maximum flux across cell surface, used in the growth 

rate of a single cell = 6 x 10 gm/cm -hr. 

K , MichaelSfe - Menten constant in the flux term of growth 
s 

rate expression = 0.02 gm/liter 

JA- f proportionality constant for rate of release of mass 
c 

by cell = 1 hr 

fraction of substrate in mass taten into the cell — 0.75 

C°. substrate concentration in entering stream =2.5 gm/cm 
s’ 

3 

p, mass density of cell = 1.01 gm/cm . 
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NOME N CLATURE 

C Biomass concentration, grams/litre 

C i Coefficient of ith approximating function in trial 

solution 

C sk Concentre tion of kth substrate (C g if single limiting 
substrate), grams/litre. 

C° k Concentration of kth substrate in feed stream, grams/litre 

f(t) Function of time in cell growth rate expression 

g*(m) Density of cell mass distribution in a symptotic growth 

-1 

gram . 

g(m) Initial cell mass distribution cells/litre-gram 

h' (m) Density of division mass distribution gram 

K g Michaelis-Menten kinetics saturation constant, gram/li 4 ‘” 

L n nth Daguerre polynomial 

1 Length of rod-like cell, cm. 

m,m* Cell mass, gram 

m c Mean division mass of a cell, gram 

N Population density, cells/litre 

Number of Apnroximating functions in trial solution. 

-1 

p(m,m T ) Density of daughter cell mass distribution, gram 

R Radius of coll (rod), cm. 

R(m) Residual of the differential equation 

r(m,C sk )C e11 growth rate grams/hr. in environment C gk 

r' Rate of cellular mass uptake grams/ljr. 

r n Rate of cellular mass release , grams/hr. 

2 

S Surface area of a cell, cm 

Time, hr. 


t 



W(m,t) Number density of cells, c ell s/gm- litre; 

w(x,t) Transformed number density cells/gm-litre 

x,y Variables used for mathematical manipulations, defined 
where used in text. 

GREEK LETTERS 




p 


r v 

r 

c 


fraction of ith component in mass taken into cell {& 
if a single limiting substrate - ) 

_ -J 

Specific probability of fission hr 

fraction of component i in mass released by cell. 

Dirac delta function 

Measure of spread in division mass distribution, gram 


@(m,C J Specific de^th probability hr~^ 


9 

x 

A- 

A 


1 

e 


m 


n 


av/b 


Holding time, hr. 

Constant in weighting function expression 

-1 

Specific rate of biomass increase hr 
ith moment of the distribution density 
Specific mass release rate hr 
Specific rate of population increase 

rz 

Mass density of cell grnm/cnr 
Cell age, hr. 

2 

Mass flux across cell surface gram/cm -hr. 

2 

Maximum mass flux across cell surface gram/ cm -hr. 
nth approximating function in trial solution 
nth weighting function 

Over a quantity means it is evaluated in the steads’- 
Over a quantity means it is a vector 
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